
!--------------------------------------------------------------------
! linearized ozone chemistry LINOZ
! from Hsu and Prather, JGR, 2008
!
! written by Jean-Francois Lamarque (September 2008)
! modified by
!     24 Oct 2008 -- Francis Vitt
!      9 Dec 2008 -- Philip Cameron-Smith, LLNL, -- added ltrop
!--------------------------------------------------------------------
module lin_strat_chem

  use shr_kind_mod, only : r8 => shr_kind_r8
  use ppgrid,       only : begchunk, endchunk
  use physics_types,only : physics_state
  use cam_logfile,  only : iulog
  use abortutils,   only : endrun
  use spmd_utils,   only : masterproc
  !
  implicit none
  !
  private  ! all unless made public

  save
  !
  ! define public components of module
  !
  public :: lin_strat_chem_inti, lin_strat_chem_solve
  public :: do_lin_strat_chem

  integer :: index_o3
  logical :: do_lin_strat_chem


contains

!--------------------------------------------------------------------
!--------------------------------------------------------------------
  subroutine lin_strat_chem_inti(phys_state)
    !
    ! initialize linearized stratospheric chemistry by reading
    ! input parameters from netcdf file and interpolate to
    ! present model grid
    !
    use linoz_data,   only : linoz_data_init, has_linoz_data
    use ppgrid,       only : pver
    use mo_chem_utls, only : get_spc_ndx
    use cam_history,  only : addfld, phys_decomp, add_default
    implicit none


    type(physics_state), intent(in) :: phys_state(begchunk:endchunk)


    if (.not.has_linoz_data) return

    !
    ! find index of ozone
    !
    index_o3 = get_spc_ndx('O3')
    do_lin_strat_chem = has_linoz_data
    if ( index_o3 <= 0 ) then
       write(iulog,*) ' No ozone in the chemical mechanism, skipping lin_strat_chem'
       do_lin_strat_chem = .false.
       return
    end if

    ! check for synoz

    if( get_spc_ndx( 'SYNOZ' ) > 0 .and. has_linoz_data) then
       call endrun('lin_strat_chem_inti: cannot have both synoz and linoz')
    endif

    ! initialize the linoz data

    call linoz_data_init()

    ! define additional output

    call addfld( 'LINOZ_DO3'    , '/s'     , pver, 'A', 'ozone vmr tendency by linearized ozone chemistry'  , phys_decomp )
    call addfld( 'LINOZ_DO3_PSC', '/s'     , pver, 'A', 'ozone vmr loss by PSCs using Carille et al. (1990)', phys_decomp )
    call addfld( 'LINOZ_SSO3'   , 'kg'     , pver, 'A', 'steady state ozone in LINOZ'                       , phys_decomp )
    call addfld( 'LINOZ_O3COL'  , 'DU'     , pver, 'A', 'ozone column above'                                , phys_decomp )
    call addfld( 'LINOZ_O3CLIM' , 'mol/mol', pver, 'A', 'climatology of ozone in LINOZ'                     , phys_decomp )
    call addfld( 'LINOZ_SZA'    , 'degrees',    1, 'A', 'solar zenith angle in LINOZ'                       , phys_decomp )

    call add_default( 'LINOZ_DO3'    , 1, ' ' )
    call add_default( 'LINOZ_DO3_PSC', 1, ' ' )
    call add_default( 'LINOZ_SSO3'   , 1, ' ' )
    call add_default( 'LINOZ_O3COL'  , 1, ' ' )
    call add_default( 'LINOZ_O3CLIM' , 1, ' ' )
    call add_default( 'LINOZ_SZA'    , 1, ' ' )

    return
  end subroutine lin_strat_chem_inti


!--------------------------------------------------------------------
!--------------------------------------------------------------------
  subroutine lin_strat_chem_solve( ncol, lchnk, o3_vmr, o3col, temp, sza, pmid, delta_t, rlats, ltrop )
 
    use chlorine_loading_data, only: chlorine_loading

    !
    ! this subroutine updates the ozone mixing ratio in the stratosphere
    ! using linearized chemistry 
    !

    use ppgrid,        only : pcols, pver
    use physconst,     only : pi, &
                              grav => gravit, &
                              mw_air => mwdry
    use cam_history,   only : outfld
    use linoz_data,    only : fields, o3_clim_ndx,t_clim_ndx,o3col_clim_ndx,PmL_clim_ndx,dPmL_dO3_ndx,&
                                      dPmL_dT_ndx,dPmL_dO3col_ndx,cariolle_pscs_ndx
    !
    ! dummy arguments
    !
    integer,  intent(in)                           :: ncol                ! number of columns in chunk
    integer,  intent(in)                           :: lchnk               ! chunk index
    real(r8), intent(inout), dimension(ncol ,pver) :: o3_vmr              ! ozone volume mixing ratio
    real(r8), intent(in)   , dimension(ncol ,pver) :: o3col               ! ozone column above box (mol/cm^2)
    real(r8), intent(in)   , dimension(pcols,pver) :: temp                ! temperature (K)
    real(r8), intent(in)   , dimension(ncol )      :: sza                 ! local solar zenith angle
    real(r8), intent(in)   , dimension(pcols,pver) :: pmid                ! midpoint pressure (Pa)
    real(r8), intent(in)                           :: delta_t             ! timestep size (secs)
    real(r8), intent(in)                           :: rlats(ncol)         ! column latitudes (radians)
    integer,  intent(in)   , dimension(pcols)      :: ltrop               ! chunk index
    !
    ! local
    !
    integer :: i,k,n !,index_lat,index_month
    real(r8) :: o3col_du,delta_temp,delta_o3col,o3_old,o3_new,delta_o3
    real(r8) :: max_sza, psc_loss
    real(r8) :: o3_clim
    real(r8), dimension(ncol) :: lats
    real(r8), dimension(ncol,pver) :: do3_linoz,do3_linoz_psc,ss_o3,o3col_du_diag,o3clim_linoz_diag

    real(r8), dimension(:,:), pointer :: linoz_o3_clim
    real(r8), dimension(:,:), pointer :: linoz_t_clim
    real(r8), dimension(:,:), pointer :: linoz_o3col_clim
    real(r8), dimension(:,:), pointer :: linoz_PmL_clim
    real(r8), dimension(:,:), pointer :: linoz_dPmL_dO3
    real(r8), dimension(:,:), pointer :: linoz_dPmL_dT
    real(r8), dimension(:,:), pointer :: linoz_dPmL_dO3col
    real(r8), dimension(:,:), pointer :: linoz_cariolle_psc

    !
    ! parameters
    !
    real(r8), parameter :: convert_to_du = 1._r8/(2.687e16_r8)      ! convert ozone column from mol/cm^2 to DU
    real(r8), parameter :: degrees_to_radians = pi/180._r8          ! conversion factors
    real(r8), parameter :: radians_to_degrees = 180._r8/pi
    real(r8), parameter :: temp_activation_cariolle = 193._r8       ! O3 loss freq when T below temp_activation_cariolle (K)
    real(r8), parameter :: chlorine_loading_1987    = 2.5977_r8     ! EESC value (ppbv)
    real(r8), parameter :: chlorine_loading_bgnd    = 0.0000_r8     ! EESC value (ppbv) for background conditions
    real(r8), parameter :: pressure_threshold       = 210.e+2_r8    ! {PJC} for diagnostics only

    !
    ! skip if no ozone field available
    !
    if ( .not. do_lin_strat_chem ) return

    ! 
    ! associate the field pointers
    !
    linoz_o3_clim      => fields(o3_clim_ndx)      %data(:,:,lchnk )
    linoz_t_clim       => fields(t_clim_ndx)       %data(:,:,lchnk )
    linoz_o3col_clim   => fields(o3col_clim_ndx)   %data(:,:,lchnk )
    linoz_PmL_clim     => fields(PmL_clim_ndx)     %data(:,:,lchnk )
    linoz_dPmL_dO3     => fields(dPmL_dO3_ndx)     %data(:,:,lchnk )
    linoz_dPmL_dT      => fields(dPmL_dT_ndx)      %data(:,:,lchnk )
    linoz_dPmL_dO3col  => fields(dPmL_dO3col_ndx)  %data(:,:,lchnk )
    linoz_cariolle_psc => fields(cariolle_pscs_ndx)%data(:,:,lchnk )

    !
    ! initialize output arrays
    !
    do3_linoz         = 0._r8
    do3_linoz_psc     = 0._r8
    o3col_du_diag     = 0._r8
    o3clim_linoz_diag = 0._r8
    ss_o3             = 0._r8
    !
    ! convert lats from radians to degrees
    !
    lats = rlats * radians_to_degrees

    LOOP_COL: do i=1,ncol
       LOOP_LEV: do k=1,ltrop(i)
          !
          ! climatological ozone
          !
          o3_clim = linoz_o3_clim(i,k)
          !
          ! skip if not in the stratosphere
          !
          if ( pmid(i,k) > pressure_threshold ) THEN   ! PJC diagnostic
             WRITE(iulog,*)'LINOZ WARNING: Exceeded PRESSURE threshold (i,k,p_threshold,pmid,o3)=',&
               i,k,nint(pressure_threshold/100.),'mb',nint(pmid(i,k)/100.),'mb',nint(o3_vmr(i,k)*1e9),'ppb'   !PJC
!             cycle LOOP_LEV
          endif
          !
          ! diagnostic for output
          !
          o3clim_linoz_diag(i,k) = o3_clim
          !
          ! old ozone mixing ratio
          !
          o3_old = o3_vmr(i,k)
          !
          ! convert o3col from mol/cm2
          !
          o3col_du = o3col(i,k) * convert_to_du
          o3col_du_diag(i,k) = o3col_du
          !
          ! compute differences from climatology
          !
          delta_temp  = temp(i,k) - linoz_t_clim    (i,k)
          delta_o3col = o3col_du  - linoz_o3col_clim(i,k)


          !
          ! steady state ozone
          !
          ss_o3(i,k) = o3_clim - (               linoz_PmL_clim   (i,k)   &
                                 + delta_o3col * linoz_dPmL_dO3col(i,k)   &
                                 + delta_temp  * linoz_dPmL_dT    (i,k)   &
                                             ) / linoz_dPmL_dO3   (i,k)


          !
          ! ozone change
          !
          delta_o3 = (ss_o3(i,k)-o3_old) * (1._r8 - exp(linoz_dPmL_dO3(i,k)*delta_t))
          !
          ! define new ozone mixing ratio
          !
          o3_new = o3_old + delta_o3
          !
          ! output diagnostic
          !
          do3_linoz(i,k) = delta_o3/delta_t
          !
          ! PSC activation (follows Cariolle et al 1990.)
          !
          ! use only if abs(latitude) > 40.
          !
          if ( abs(lats(i)) > 40. ) then   
             if ( (chlorine_loading-chlorine_loading_bgnd) > 0. ) then
                if ( temp(i,k) <= temp_activation_cariolle ) then
                   !
                   ! define maximum SZA for PSC loss (= tangent height at sunset)
                   !
                   max_sza = (90._r8 + sqrt( max( 16._r8*log10(100000._r8/pmid(i,k)),0._r8)))
#ifdef DEBUG
                   write(iulog,*)sza(i),max_sza
#endif
                   if ( (sza(i)*radians_to_degrees) <= max_sza ) then

                      psc_loss = exp(-linoz_cariolle_psc(i,k) &
                           * (chlorine_loading/chlorine_loading_1987)**2 &
                           * delta_t )

                      o3_new = o3_old * psc_loss
                      !
                      ! output diagnostic
                      !
                      do3_linoz_psc(i,k) = (o3_new-o3_old)/delta_t
                      !
                   end if
                end if
             end if
          end if
          !
          ! update ozone vmr
          !
          o3_vmr(i,k) = o3_new

       end do LOOP_LEV
    end do LOOP_COL
    !
    ! output
    !
    call outfld( 'LINOZ_DO3'    , do3_linoz              , ncol, lchnk )
    call outfld( 'LINOZ_DO3_PSC', do3_linoz_psc          , ncol, lchnk )
    call outfld( 'LINOZ_SSO3'   , ss_o3                  , ncol, lchnk )
    call outfld( 'LINOZ_O3COL'  , o3col_du_diag          , ncol, lchnk )
    call outfld( 'LINOZ_O3CLIM' , o3clim_linoz_diag      , ncol, lchnk )
    call outfld( 'LINOZ_SZA'    ,(sza*radians_to_degrees), ncol, lchnk )

    return
  end subroutine lin_strat_chem_solve

end module lin_strat_chem
